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ABSTRACT 

A  methodology  for  choosing  a  hybrid  cr-p  (sigma-pressure)  vertical  coordinate  of  the  Simmons-Strtifing 
form  for  a  global  model  is  presented.  The  method  focuses  on  properties  of  the  vertical  derivative  of  the 
terrain-following  coefficient,  which  affect  the  smoothness  and  shape  of  layer  thickness  profiles  and  deter¬ 
mines  the  coordinate’s  nronotonicity  over  variable  terrain.  The  method  is  applied  to  characterize  and 
interrelate  existing  hybrid  coordinate  choices  in  NWP  and  climate  models,  then  to  design  new  coordinates 
with  specific  properties.  Offline  tests  indicate  that  the  new  coordinates  reduce  stratospheric  errors  in  models 
due  to  vertical  truncation  effects  in  the  computation  of  the  pressure  gradient  force  over  steep  terrain.  When 
implemented  in  a  global  model,  the  new  coordinates  significantly  reduce  vorticity  and  divergence  errors  at 
all  altitudes  in  idealized  simulations.  In  forecasting  experiments  with  a  global  model,  the  new  coordinates 
slightly  reduce  the  stability  of  the  semi-implicit  time  scheme.  Resetting  the  reference  pressure  in  the  scheme 
to  —800  hPa  solves  the  problem  for  every  coordinate  except  the  Sangster-Arakawa-Lamb  hybrid,  which 
remains  intrinsically  less  stable  than  the  others.  Impacts  of  different  coordinates  on  forecast  skill  are  neutral 
or  weakly  positive,  with  the  new  hybrid  coordinates  yielding  slight  improvements  relative  to  earlier  hybrid 
choices.  This  essentially  neutral  impact  indirectly  endorses  the  wide  variety  of  hybrid  coordinate  choices 
currently  used  in  NWP  and  climate  models,  with  the  proviso  that  these  tests  do  not  address  the  impact  over 
longer  time  scales  or  on  data  assimilation. 


1.  Introduction 

Many  global  models  solve  the  hydrostatic  primitive 
equations  on  the  sphere  using  a  generalized  vertical 
coordinate,  based  on  equation  sets  first  derived  in  con¬ 
tinuous  form  by  Kasahara  (1974)  and  in  discretized 
form  by  Simmons  and  Burridge  (1981).  Given  pressure 
p  as  the  physical  altitude  variable,  this  coordinate  can 
take  the  functional  form: 

t l  =  h(p,ps),  (1) 

where  ps  is  the  surface  pressure  and  h(p,ps )  is  a  mono¬ 
tonic  function  of  p,  such  that  there  is  a  one-to-one  map¬ 
ping  between  p  and  tj  for  a  given  ps.  The  functional 
dependence  on  ps  permits  terrain-following  coordi¬ 
nates.  The  standard  convention  is  to  define  tj  to  be 
normalized  and  dimensionless,  satisfying  the  lower  and 
upper  boundary  conditions: 
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h(ps,Ps)  =  1, 

(2) 

h(P  toV,  Ps)  =  0, 

(3) 

respectively,  where  ptop  is  the  pressure  at  the  top  model 
half-layer.  The  chain  rule  applied  to  (1)  allows  rj  terms 
in  the  equations  to  be  reexpressed  in  terms  of  p ,  ps ,  and 
ptop,  which  has  the  practical  advantage  that  rj  need  not 
be  defined  explicitly. 

The  most  straightforward  choice  for  rj  is  the  terrain¬ 
following  cr  coordinate,  obtained  by  setting 

V )  =  h(p,  ps)  =  '  _/  top  =  o-.  (4) 

Ps  P  top 

This  is  a  minor  generalization  of  the  original  Phillips 
(1957)  cr  coordinate  pip s  for  models  with  a  nonvanish¬ 
ing  Ptop  (Mintz  1965). 

However,  as  models  extend  into  the  middle  atmo¬ 
sphere,  there  are  numerical  and  practical  advantages  in 
choosing  a  function  h(p,ps )  that  smoothly  transitions  rj 
from  terrain-following  levels  near  the  ground  to  iso- 
baric  surfaces  in  the  upper  troposphere  or  lower  strato¬ 
sphere  (Fels  et  al.  1980;  Simmons  and  Burridge  1981; 
Simmons  and  Striding  1983;  Simmons  et  al.  1989).  Sim- 
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mons  and  Striding  (1981)  and  Simmons  and  Burridge 
(1981)  tested  different  functional  forms  for  h(p ,  ps) 
that  yielded  so-called  hybrid  cr-p  coordinates.  The  final 
function  considered  by  Simmons  and  Striding  (1981) 
took  the  implicit  form 


p(r ),  Ps)  =  Mv)p0  +  B(ff)ps,  (5) 

where  p0  is  some  nominal  sea  level  pressure,  typically 
—1000  hPa,  and  fj  is  the  corresponding  vertical  profile 
of  Tj  values  for  ps  =  p0.  Later,  Simmons  and  Striding 
(1983)  replaced  A(f))p0  with  the  simpler  form  A(f)). 
Generalizing  to  a  nonzero  plop  yields 

PiV’Ps)  =  A(v)  +  B(v)(Ps  ~  Ptop ),  (6) 

to  satisfy  the  upper-level  boundary  condition  in  (3).  A 
corresponding  explicit  expression  for  t)  that  satisfies  the 
boundary  conditions  (2)  and  (3)  is 


P  ~  Ptop 
PS  ~  Ptop 


My)  ~  Ptop 
PS  ~  Ptop 


+  B(fj), 


(7) 


such  that  /;lop  <  A(f))  <  p,  0  <  B( fj)  <  1,  and  dBldr\  >  0. 

This  coordinate  is  defined  by  vertical  profiles  of  the 
two  coefficients,  A  and  B.  which  control,  respectively, 
its  isobaric  and  terrain-following  properties.  For  B  J=  0, 
the  coordinate  has  terrain-following  properties,  from 
(6),  of  the  following  form: 


B(rj) 


P  ~  Mv) 
PS  ~  Pt op  ' 


(8) 


By  analogy  to  (4),  (8)  implies  a  cr-like  coordinate  on  the 
reference  level  f)  in  which  the  effective  model  top  is  not 
pto p,  but  A{  fj).  Since  a  surfaces  of  the  form  in  (4)  flatten 
out  as  p  — >  ptop  (see  Fig.  1  of  Hogan  and  Rosmond 
1991)  and  A(tj)  >  ptop,  (8)  yields  a  surface  whose  re¬ 
sponse  to  terrain  has  been  reduced,  or  flattened,  rela¬ 
tive  to  pure  <j  (4)  when  A(fj)  >  ptop. 

While  more  elaborate  functional  forms  for  /;(/?,  ps) 
can  be  chosen  (e.g.,  Simmons  and  Burridge  1981;  Ar- 
akawa  and  Konor  1996),  the  Simmons  and  Strafing 
(1983)  form  in  (6)  has  proven  sufficiently  flexible  that  it 
is  now  built  into  the  dynamical  cores  of  many  numerical 
weather  prediction  (NWP)  and  climate  models.  Fur¬ 
thermore  (6)  can  be  used,  with  minor  modifications,  in 
nonhydrostatic  dynamical  cores  (Laprise  1992;  Wood 
and  Staniforth  2003),  and  can  be  extended  to  accom¬ 
modate  more  sophisticated  hybrid  formulations,  such 


as  isentropic  levels  0  through  a  hybrid  cr-0-p  coordinate 
of  a  similar  functional  form  in  which  h  is  now  also  a 
function  of  temperature  (e.g.,  Zhu  et  al.  1992;  Webster 
et  al.  1999). 

Yet,  despite  its  widespread  use  in  NWP  and  climate 
models,  there  is  no  specific  guidance  in  the  literature  on 
how  best  to  set  the  governing  coefficients  A(p)  and 
B( fj)  to  provide  an  appropriate  hybrid  coordinate  for  a 
given  modeling  application.  Simmons  and  Burridge 
(1981),  who  presented  subsets  of  the  results  in  Simmons 
and  Striifing  (1981),  did  not  consider  hybrid  coordi¬ 
nates  of  the  form  in  (6)  at  all  (apart  from  a  brief  foot¬ 
note  added  in  proof).  While  presenting  a  series  of  im¬ 
portant  model  tests  relevant  to  the  choice  of  a  hybrid 
coordinate,  Simmons  and  Burridge  (1981)  explicitly 
made  no  specific  recommendations  about  which  hybrid 
coordinate  to  use,  noting  choices  were  dependent  on 
specific  modeling  application  and  spatial  resolution. 
The  global  modeling  study  of  Simmons  and  Striifing 
(1983)  compared  several  different  hybrid-coordinate 
and  cr-coordinate  forecasts,  finding  small  forecasting 
improvements  with  hybrid  coordinates  and  limited  dif¬ 
ferences  in  model  performance  among  the  various  hy¬ 
brid  coordinates  they  tested.  Again,  no  specific  guid¬ 
ance  on  an  optimal  choice  for  A{ r))  and  /1(f))  was  pro¬ 
vided,  apart  from  one  hybrid  coordinate  which 
rendered  their  model  more  unstable  numerically  and 
was  discarded.  Due  perhaps  to  the  small  effects  on  fore¬ 
cast  skill  noted  by  Simmons  and  Striifing  (1983),  there 
have  been  no  subsequent  tests  of  hybrid  cr-p  coordi¬ 
nates  reported  in  the  literature. 

The  main  purpose  of  this  paper  is  to  set  forth  a 
straightforward  methodology  for  designing  and  testing 
hybrid  cr-p  vertical  coordinates  of  the  Simmons  and 
Striifing  (1983)  form  for  a  global  model.  In  section  2, 
we  show  that  the  coordinate’s  relevant  properties,  such 
as  smoothness  and  shape  preservation  of  layer  thick¬ 
ness  profiles  and  monotonicity  over  high  terrain, 
are  encapsulated  in  the  vertical  gradient  of  the  terrain¬ 
following  coefficient,  dBldr).  We  apply  these  dBldr) 
diagnostics  in  section  3  to  introduce  and  characterize  a 
series  of  hybrid  vertical  coordinates  used  previously  in 
models,  but  now  all  expressed  in  the  Simmons  and 
Striifing  (1983)  form  in  (6).  In  section  4  we  use  the  same 
diagnostics  to  design  some  new  hybrid  coordinates 
with  specific  properties.  Section  5  studies  the  errors 
in  the  pressure  gradient  force  produced  by  these  coor¬ 
dinates  over  steep  terrain  due  to  truncation  effects  of 
the  global  model’s  vertically  discretized  finite-dif¬ 
ference  formulation.  Section  6  implements  and  tests  hy¬ 
brid  coordinates  in  an  advanced-level  physics  high- 
altitude  (ALPHA)  prototype  of  the  Navy  Operational 
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Global  Atmospheric  Prediction  System  (NOGAPS),  to 
investigate  their  effects  on  model  errors,  stability,  and 
forecast  skill.  Our  findings,  along  with  coordinate 
choices  used  in  other  climate  and  NWP  models,  are 
summarized  in  section  7. 

2.  Formulation  and  diagnostics 

a.  Discretization 

In  cr-coordinate  models  it  is  common  to  tag  model 
levels  by  their  invariant  cr  values.  This  fails  for  hybrid 
coordinates:  from  (7),  any  isobaric  component  to  the 


coefficient  [A( fj)  >  ptop]  causes  tj  values  to  vary  with 
surface  pressure.  Thus,  we  tagged  model  levels  initially 
in  section  1  using  the  profile  rj  evaluated  at  a  constant 
nominal  sea  level  pressure  p0,  such  that  A(fj)  and  B( fj) 
are  the  invariant  values  of  the  isobaric  and  terrain¬ 
following  coefficients,  respectively,  on  a  given  model 
level  fj. 

In  discretized  form,  full  model  levels  (layers)  are 
tagged  from  top  to  bottom  with  indices  k  =  1,  .  . . ,  L, 
where  L  is  the  total  number  of  layers.  Pressures  on  the 
L  +  1  half  (interface)  levels  are  specified  using  the 
hybrid  cr-p  function  (6): 


P(Vk+v2,Ps)  =  A(Vk+ 1/2)  +  B(r)k+1/2)(ps  ~  /hop),  now  equivalently 

Pk+l/2  =  Ak+ 1/2  +  Bk+m  (ps  —  Ptop)-  (9) 


To  reproduce  ptov  and  ps  at  the  top  and  bottom  half 
levels,  respectively,  the  boundary  conditions  for  (9), 
analogous  to  (2)  and  (3),  are 

A(0)  =  AV2  =  ptop  B(  0)  =  B1/2  =  0, 

A(l)  =  Al+1/2  =  ptop  5(1)  =  BL+lf2  =  1.  (10) 


In  the  limit  of  an  infinitesimal  surface  pressure  change, 
(15)  can  be  reexpressed  as 


MPk 

dPs 


\Bk. 


(16) 


Layer  thicknesses  are 

&Pk  =  Pk+1/2  ~  Pk- 1/2  =  AA*  +  A Bk{ps  -  ptop), 

(11) 

where 


b.  Shape  preservation  of  layer  thickness  profiles 

Vertical  profiles  of  pressure  thicknesses  A pk  are  care¬ 
fully  chosen  in  global  models.  Figure  1  plots  represen¬ 
tative  profiles  of  A pk  and  the  corresponding  pressure- 
height  thicknesses  A Zk,  the  latter  computed  from  pres¬ 
sure  altitudes, 


—  Ak+1/2  Ak_y 2, 

A Bk  =  Bk+lf2  ~  Bk_1/2.  (12) 

From  (7),  (12)  also  defines  the  model’s  t j  layer  thick¬ 
nesses: 

A Ak 

=  ~  +  A Bk.  (13) 

PS  P  top 

The  fj  levels,  evaluated  at  a  nominal  sea  level  surface 
pressure  p0  ~  1000  hPa,  yield  corresponding  half-level 
pressures  pk+u 2,  such  that  pL+y2  =  p0.  Consider  now  a 
surface  pressure  ps  that  differs  from  this  nominal  sea 
level  value  pL+1/2.  Using  (9)  we  can  express  the  new 
pressure  levels  pk+ V2  and  layer  pressure  thicknesses 
A pk  as  departures  from  the  nominal  profiles  pk+ 1/2  and 

Apfc: 

Pk+1/2  =  Pk+1/2  +  Bk+1/2(ps  —  pL+ 1/2),  (14) 

fyk  =  ^Pk  +  ^Bk(ps  -  pL+V2).  (15) 


Zk+ia  =  (17) 

\P  k+1/2  ) 


using  Pq  =  1000  hPa  and  H  =  1  km.  The  different 
curves  compare  the  operational  NOGAPS  L30  and  re¬ 
search  NOGAPS-ALPHA  L60  and  L84  thicknesses 
with  some  other  operational  configurations  used  in  the 
European  Centre  for  Medium-Range  Weather  Fore¬ 
casts  (ECMWF)  Integrated  Forecast  System  (IFS)  and 
National  Centers  for  Environmental  Prediction 
(NCEP)  Global  Forecast  System  (GFS).  In  addition  to 
concentrating  layers  in  preferred  regions  such  as  the 
boundary  layer,  these  profiles  vary  thicknesses 
smoothly  across  adjacent  layers  in  to  order  to  minimize 
introduction  of  errors  through  the  model’s  vertical  fi¬ 
nite-difference  calculations. 

Equation  (16)  determines  how  these  thickness  pro¬ 
files  respond  over  terrain.  We  assume  choices  for  B{ fj) 
that  vary  B  values  sufficiently  gradually  in  comparison 
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(a)  Pressure  Layer  Thicknesses 


(b)  Pressure  Height  Layer  Thicknesses 


Oh 

a 


Fig.  1.  Plot  of  layer  thicknesses  (a)  A pk  and  (b)  A Zk  vs  pk  and  Zk,  respectively.  Results  are  shown  for  two  recent 
ECMWF  IFS  configurations  (L91,  L60),  the  operational  NOGAPS  L30  levels,  two  sets  of  NOGAPS-ALPHA 
levels  (L60,  L84),  and  a  recent  set  of  operational  NCEP  GFS  L64  levels. 


to  model  layer  thicknesses  that  the  linear  approxima¬ 
tion 

(dB\ 

(18) 

is  accurate,  whereupon  (15)  can  be  reexpressed  as 
(dB\ 

fyk  ~  ^Pk  +  J  Arjk(Ps  ~  pL+in)-  (19) 

Since  Ar)fc  =  A pkl(pL+v2  ~  Ptop).  we  can  reexpress 
(19)  as 


pressures  that  yield  vanishing  thicknesses  A pk,  which 
from  (22)  occur  at  any  level  f)  where  s(r]k,ps)  =  0. 
From  (21),  this  latter  criterion  can  be  expressed  equiva¬ 
lently  as 

dB/dr)  -  1  =  (ps)m -  ptop 
dB/drj  pL+ ^  ~  ptop  ' 

For  any  5(r )),  (ps)min  in  (23)  is  the  minimum  surface 
pressure  that  this  coordinate  can  tolerate  before  be¬ 
coming  nonmonotonic. 

c.  Setting  A(r))  and  5(f)) 


Apk  =  ^Pk 


dB 
1  +  I  “TT 


&Ps 


dv  /  k\pL  + 1/2  Pto 


(20) 


where  A ps  =  ps  —  pL+y2  is  the  change  in  surface  pres¬ 
sure  from  the  nominal  sea  level  value  pL+y2-  The  mul¬ 
tiplicative  term  in  square  parentheses  in  (20): 


S^Ps)  =  1  +  [dtJ 


dB\ 


k  \P L  +  l/2  P to 


(21) 


is  a  “shape  factor”  that  modifies  the  values  and  shape  of 
the  reference  thickness  profile  A pk  over  terrain,  such 
that  the  new  profile  is 


&pk  =  kpks(T\k,ps).  (22) 


From  (22),  we  get  perfect  shape  preservation  of  pres¬ 
sure  thickness  profiles  when  ds/d f)  =  0,  or  equivalently 
from  (21)  when  dBldr)  is  a  constant. 

The  definition  in  (1)  requires  tj  to  be  a  monotonic 
function  of  p  for  all  possible  surface  pressures  ps  in 
the  model.  We  can  check  this  by  testing  for  surface 


The  primary  importance  of  dB/dr]  in  defining  the  sa¬ 
lient  characteristics  of  hybrid  coordinates  of  the  Sim¬ 
mons  and  Strafing  (1983)  form  (6)  implies  that  its  re¬ 
quired  form  should  drive  algorithms  for  setting  the  co¬ 
efficients  Ak+m  and  Bk+y2  in  a  model.  A  flexible  and 
fairly  general  algorithm  for  assigning  Ak+V2  and  Bk+V2 
values  based  on  this  philosophy  is  described  in  the  ap¬ 
pendix. 

The  algorithm  allows  the  top  kp  model  layers  to  be 
isobaric  and  lowest  k,r  model  layers  to  be  cr-like.  The 
intervening  model  layers  k  =  kp  +  1, . . . ,  L  —  k„  take 
a  hybrid  form  defined  by  terrain-following  coefficients 
of  the  analytical  form: 


5(7)) 


/v~  Vkp+i/2\rm 
V1  "  Pkp  +  l/2/ 


b(r])rm. 


(24) 


The  A( tj)  values  follow  from  (7)  using  the  reference 
values  tj  =  t)  and  ps  =  pL+ll2. 

The  terrain-following  coefficient  gradient  is  given  by 
the  analytical  expression: 
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(e)  NOGHYB 
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Fig.  2.  T239L84  half-levels  pk+1/2  at  ~34.7°N  from  Zk+ll2  =  0-40  km  for  all  vertical  coordinates  tested  in  NOGAPS-ALPHA  (see 
Table  1):  (a)  cr,  (b)  SAL,  (c)  LG,  (d)  CMAM,  (e)  original  NOGAPS-ALPHA  hybrid  (NOGHYB),  and  (f),  (g)  two  new  NOGAPS- 
ALPHA  hybrids  (NEWHYB1,  NEWHYB2).  The  pk+ V2  contours  are  color  coded,  with  light  (dark)  blue  lines  indicating  Bk+1/2  near 
zero  (unity).  Green  curves  show  the  lowest-altitude  isobaric  level  pkp+y2  =  86.7  hPa  ( kp  =  58)  for  the  coordinates  with  kp  >  0.  T239 
topographic  elevations  are  shaded  light  brown. 


dB 

drj 


'0 

r(rj)b(7j) 


r(  n)-l 


1  - 


Vkp  +  1/2. 

.  (1  ~  Bkp+ 1/2) 


+  b{rj) 


*(l) 


r(n) 


k  =  l  ,...,kp, 
k  =  kp  +  1, . . . ,  L  -  k„ 
k  =  L  +  1  —  ka, . . . ,  L. 


(25) 


Thus,  the  form  of  this  hybrid  coordinate,  defined  by 
(25),  is  dictated  solely  by  the  choices  for  kp,  k (r,  and  the 
power-law  index  profile  r( f)),  which  are  varied  until  a 
coordinate  with  the  requisite  dBldr)  properties  is  ob¬ 
tained. 

3.  Examples 

We  use  the  algorithm  described  in  section  2c  to  in¬ 
troduce  some  basic  hybrid  coordinates,  which  we  then 


characterize  using  the  rfB/t/rj-based  diagnostic  tools  of 
section  2b. 

a.  cr  coordinate 

The  cr  coordinate  [(4)]  is  recovered  by  setting  kp  =  0 
and  k,r  =  L,  or  equivalently  by  setting  kp  =  k,r  =  0  and 
r(Vk r+1/2)  =  rk+m  =  1-  Resulting  half-levels  around  a 
34.7°N  latitude  circle  are  plotted  in  Fig.  2a  for  a 
T239L84  NOGAPS-ALPHA  model  resolution  (see 
Fig.  1). 


pressure  (hPa) 
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dB/dx\ 

Fig.  3.  Profiles  of  dB!dr\  for  the  vertical  coordinates  discussed 
in  section  3  (see  Table  1). 

Our  diagnostic  tools  of  section  2b  illustrate  some  of 
the  advantages  of  this  coordinate.  Since  A(t))  =  ptop 
and  B{i j)  =  fj,  then  dBldf]  =  1  and  ds/d fj  =  0  for  all  fj. 
Thus,  the  a  coordinate  perfectly  preserves  the  shapes  of 
the  model  thickness  profiles  A pk  in  Fig.  la.  Furthermore, 
the  monotonicity  criterion  in  (23)  yields  (p_s)min  =  ptop. 
Thus,  the  coordinate  is  always  monotonic  so  long  as 
mountain  peaks  do  not  breach  the  top  of  the  model. 

The  shape-preserving  dBldf]  =  1  profile  of  the  cr  co¬ 
ordinate  is  plotted  in  Fig.  3. 

b.  Sangster-Arakawa-Lamb  coordinate 

This  hybrid  coordinate  is  used  in  a  number  of  climate 
models  (e.g.,  Cariollc  and  Deque  1986;  Eluszkiewicz  et 
al.  2000;  Collins  et  al.  2004),  and  is  perhaps  the  most 
straightforward  to  implement  and  visualize.  In  the  ter¬ 
rain-following  region  k  >  kp ,  it  has  the  analytical  form 
iP  ~  Pkp+i/2)/(Ps  ~  Pkp+ 1/2)  (Arakawa  and  Lamb  1977). 
Thus,  it  consists  of  a  cr  coordinate  in  (4)  formulated 
between  the  surface  and  pk  +V2 ,  with  isobaric  model 
layers  stacked  on  top  between  pk  +1/2  and  ptop.  While 
Haltiner  and  Williams  (1980)  and  Simmons  and  Burr- 
idge  (1981)  attribute  this  coordinate  to  Arakawa  and 
Lamb  (1977),  who  derived  it  for  use  in  the  University  of 
California,  Los  Angeles  (UCLA)  general  circulation 
model  (Schlesinger  and  Mintz  1979),  it  was  derived  ear¬ 
lier  by  Sangster  (1960),  whose  ideas  were  used  to  build 
this  coordinate  into  the  Geophysical  Fluid  Dynamics 
Laboratory  (GFDL)  SKYFII  model  (Fels  et  al.  1980). 
Thus,  here  we  refer  to  it  as  the  Sangster-Arakawa- 
Lamb  (SAL)  coordinate.  Simmons  and  Strufing  (1981) 
referred  to  it  as  the  “traditional”  hybrid  coordinate. 

For  models  using  the  Simmons  and  Strufing  (1983) 


(a)  160°E:  ps=  lOOOhPa  (b)  0°E:  ps=  859hPa 


Apt  (hPa)  Apt  (hPa) 


(c)  1 10°W:  p,=  758hPa  (d)  90°E:  p,,=  469hPa 


0  20  40  60  0  20  40  60 

A Pi  (hPa)  Apk  (hPa) 


Fig.  4.  Solid  curve  shows  response  of  L84  layer  pressure  thick¬ 
nesses  A pk  vs  pk  for  the  SAL  coordinate  at  four  selected  longi¬ 
tudes  in  Fig.  2b  with  progressively  higher  terrain:  (a)  160°E,  (b) 
0°E,  (c)  110°W,  and  (d)  90°E.  Surface  pressures  ps,  given  in  plot 
title,  are  calculated  as  p0exp(—h'IH)  where  h'  is  topographic  el¬ 
evation,  H  =  7  km,  and  p0  =  1000  hPa.  The  dotted  curve  in  each 
panel  shows  thicknesses  A pk  for  ps  =  p0.  Gray  line  shows  the  lowest 
pure  pressure  interface  level  at  pkp+y2  86.7  hPa  (kp  =  58). 

form  in  (6),  it  is  obtained  by  setting  rk+1/2  =  1  in  (24). 
The  terrain-following  layers  k  >  kp  then  assume  SAL 
forms  for  any  k(T  in  the  range  0  <  k,T  <  L  —  kp  (see  the 
appendix)  of  Bk+ll2  =  bk+V2  and  Ak+V2  =  pkp+V2  + 
bk+miPtop  ~  Pkp+ 1/2)-  Figure  2b  plots  these  SAL  half 
levels  using  the  L84  reference  levels  in  Fig.  la  with  k  = 
58,  which  yields  a  lowest  isobaric  half-level  of  pk  +1/2  = 
86.7  hPa  (shown  in  green).  From  (25),  the  terrain¬ 
following  coefficient  derivative 

dB  f 0  k=l,...,kp, 

\(1  —  Pkp+1/2)  k  =  kp  +  1, . . . ,  L, 

(26) 

which  is  plotted  as  the  dotted  curve  in  Fig.  3.  Insert¬ 
ing  (26)  into  the  monotonicity  criterion  (23)  yields 

(P  s)  m  i  n  Pkp+ 1/2’ 

The  constant  terrain-following  coefficient  gradients 
in  (26)  imply  constant  shape  factors  s{r\k,ps)  and  hence 
shape  preservation  of  pressure  thicknesses  A pk,  with  one 
important  exception:  a  sharp  discontinuity  in  dBldr) 
at  the  kp  +  1/2  interface  (see  Fig.  3).  Profiles  of  A pk  in 
Fig.  4,  taken  at  selected  longitudes  in  Fig.  2b  with  pro¬ 
gressively  higher  terrain,  reveal  the  resulting  effect  on 
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layer  thicknesses  in  the  L84  model.  While  the  cr-like 
nature  of  the  coordinate  at  lower  altitudes  yields 
perfect  shape  preservation,  this  dBldi )  change  across 
the  kp  +  1/2  interface  produces  a  large  A pk  discontinu¬ 
ity  over  high  terrain  (see  Fig.  4d).  The  discontinuity  is 
even  more  apparent  in  the  corresponding  pressure- 
height  thicknesses  A Zk  in  Fig.  5.  Such  discontinuities 
can  yield  large  errors  in  the  model’s  vertical  finite- 
difference  calculations,  potentially  triggering  erroneous 
gravitational  instabilities  over  high  terrain  and  greater 
computational  instability  of  the  semi-implicit  time 
scheme  (Simmons  and  Burridge  1981).  We  return  to 
the  issue  of  model  stability  in  section  6b. 

c.  Laprise-Girard  coordinate 

This  coordinate,  developed  by  Laprise  and  Girard 
(1990)  for  the  Canadian  Climate  Center  (CCC)  model 
(McFarlane  et  al.  1992),  is  recovered  by  setting  rk+ 1/2  = 
2.  While  Laprise  and  Girard  (1990)  did  not  impose  any 
isobaric  layers  ( kp  =  0),  it  is  straightforward  to  gener¬ 
alize  to  nonzero  kp.  To  avoid  a  thickness  discontinuity 
across  the  L  —  k^  +  1/2  interface  analogous  to  that  for 
the  SAL  coordinate  at  kp  +  1/2,  we  must  set  =  0 
here. 


(a)  160°E:  ps  =  lOOOhPa 


(b)  0°E:  ps=  8591iPa 


A Zk  (km)  A Zk  (km) 

Fig.  5.  As  in  Fig.  4,  but  profiling  the  SAL-coordinate  response  of 
L84  pressure-height  layer  thicknesses  A Zk  vs  Zk  (solid)  relative  to 
the  L84  profile  A Zk  at  ps  =  p0  (dotted;  see  Fig.  lb). 

From  (25),  these  choices  yield 


dB  \  0  k  =  l,...,k 

—  —  i  (27) 

dd)  (2(fj  —  fj^+1/2)(l  —  fikp+1/2)  k  =  kp  +  1, . . . ,  L, 

which  is  plotted  as  the  dashed  curve  in  Fig.  3  using  kp  =  0.  Substituting  the  latter  expression  in  (27)  into  the 
monotonicity  criterion  (23)  yields 


[(P s)mm\k+l/2  Ptop  (Pl  +  1/2  Ptop)  1  ^ 


Vkn  + 1/2 


-  0.5 


Vk  +  l/2  Vkp  +  l/l 


(28) 


which  maximizes  at  the  surface  (f)i+1/2  =  1)  to  yield 

[(Ps)min]z.  +  l/2  =  Q-$Pl  +  l/2i  (29) 

(Laprise  and  Girard  1990).  For  pL+ 1/2  =  p0  —  1000  hPa, 
this  yields  vanishing  layer  thicknesses  at  the  surface 
when  surface  pressures  over  high  terrain  drop  below 
-500  hPa. 

Figure  2c  plots  Laprise-Girard  (LG)  coordinate  sur¬ 
faces  for  the  NOGAPS-ALPHA  L84  thicknesses  in  Fig. 
la.  The  coordinate  surfaces  transition  naturally  to 
quasi-isobaric  levels  in  the  stratosphere,  explaining  why 
Laprise  and  Girard  (1990)  implemented  it  without  an 
explicit  transition  to  pure  pressure  levels  ( kp  =  0). 
However,  over  the  Himalayas,  where  our  idealized  sur¬ 
face  pressure  calculation  yields  ps  =  469  hPa,  the  thick¬ 


nesses  in  Fig.  2c  become  negative  at  the  surface,  yield¬ 
ing  subterranean  model  levels,  consistent  with  (29). 
This  coordinate  is  evidently  stable  in  climate  models 
due  to  their  low  horizontal  resolutions,  which  reduce 
peak  terrain  heights  and  increase  minimum  terrain 
pressures.  Clearly  this  coordinate  is  more  prone  to  vio¬ 
lating  monotonicity  in  models  with  higher  horizontal 
resolution. 

This  probably  explains  the  recent  transition  of  the 
higher-resolution  Canadian  Middle  Atmosphere 
Model  (CMAM),  a  CCC  climate  model  derivative, 
from  rk+V2  =  2  to  rk+1/2  =  1.5  (Polavarapu  et  al.  2005), 
which  reduces  {ps)miJpL+i/2  from  1/2  to  1/3.  The  result¬ 
ing  dBldr\  profile  is  plotted  as  the  dot-dashed  curve  in 
Fig.  3.  The  corresponding  CMAM  coordinate  surfaces 
are  plotted  in  Fig.  2d,  showing  that  the  improved  mono- 
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tonicity  of  this  modified  coordinate  comes  at  the  ex¬ 
pense  of  a  more  gradual  transition  with  altitude  from 
terrain-following  to  isobaric  levels. 

d.  Initial  NOGAPS-ALPHA  (NOGHYB) 
coordinate 

NOGAPS-ALPHA  is  a  high-altitude  prototype  of 
the  NOGAPS  NWP  system  (Allen  et  al.  2006).  One 
aspect  of  that  extension  was  to  replace  the  cr  levels  used 
operationally  with  hybrid  levels  that  transitioned  to 
pressure  surfaces  in  the  new  middle-atmosphere  layers 
(Eckermann  et  al.  2004).  To  have  as  little  impact  as 
possible  on  NOGAPS  tropospheric  NWP  performance, 
the  goal  initially  was  to  use  model  levels  from  —1000  to 
400  hPa  that  were  as  close  as  possible  to  the  operational 
model’s  pure  cr  levels,  before  transitioning  to  newer 
hybrid  cr-p  layers  at  higher  altitudes.  These  choices  are 
reflected  in  the  thickness  profiles  in  Fig.  1:  the 
NOGAPS-ALPHA  L60  and  L84  thickness  profiles 
largely  track  the  operational  L30  thicknesses  from  1000 
to  400  hPa,  departing  in  shape  thereafter. 

The  hybrid  coordinate  that  best  retains  cr-like  tropo¬ 
spheric  layers  is  the  SAL  coordinate,  but  its  major 
drawback,  highlighted  in  section  3b,  is  the  sharp  dis¬ 
continuity  in  layer  thicknesses  over  high  terrain  at  kp  + 
1/2  (Fig.  4)  and  the  potential  impact  on  model  errors 
and  stability  (Simmons  and  Strtifing  1983).  Thus,  the 
hybrid  coordinate  implemented  in  NOGAPS-ALPHA 
used  the  algorithm  outlined  in  section  2c  with  a  power 
index  of  the  following  form: 


r(rj)  =  rp  +  (^clanP<!  arctan[5,h(fj)],  (30) 

where  r,T  is  the  limiting  index  at  L  —  ka  +  1/2,  rp  is  the 
limiting  index  at  kp  +  1/2,  and  S  is  a  dimensionless 
constant  that  controls  the  smooth  transition  to  each 
limiting  value,  such  that 


Fig.  6.  NOGHYB  coordinate  shape  factors  s(rik,ps)  for  L84 
hybrid  levels  in  Fig.  2e  at  the  indicated  longitudes  with  progres¬ 
sively  higher  terrain.  Thick  gray  curve  marks  pk  +V2  =  86.7  hPa 
(kp  =  58). 


to  avoid  SAL-like  thickness  discontinuities  across  the 
kp  +  1/2  interface. 

The  L84  NOGHYB  levels  are  plotted  in  Fig.  2e.  Fig¬ 
ure  6  plots  the  L84  shape  factors  s(yk,ps)  at  various 
longitudes  in  Fig.  2e  over  progressively  higher  terrain. 
Note  the  roughly  constant  SAL-like  shape  factors  at  0.4 
s  fj  <  1  and  the  rapid  transition  to  invariant  shapes  at 
yk  +i/2-  Figure  7  plots  the  response  of  the  L84  layer 
thicknesses  A pk.  While  this  coordinate  produces  a  sec¬ 
ondary  thickness  minimum  over  very  high  terrain  (Fig. 
7d),  the  profiles  remain  smooth  and,  in  comparison  to 
the  corresponding  plots  for  the  SAL  coordinate  in  Fig. 
4,  maintain  smoothness  and  continuity  of  the  thickness 
profiles  at  all  model  levels  k. 


drG l)  =  ~  rp )  S _ 1  , 

dy  arctanS  {l  +  [Sb(fj)]2}  (1  -  ykp+m)  ’ 

for  use  in  (25).  For  kiT  #  0,  we  require  r(r  =  1  to  avoid 
thickness  discontinuities  across  the  L  —  k^  +  1/2  inter¬ 
face. 

We  eventually  settled  upon  rp  =  2,  /• =  1,  and  S  = 
10.  These  NOGHYB  coordinate  settings  are  listed  in 
Table  1  and  its  dBldy  profile  is  plotted  as  the  solid 
curve  in  Fig.  3.  For  pk  +1/2  —  100  hPa,  these  choices 
yield  a  SAL-like  coordinate  (r  «=  1)  at  1000  hPa  S 
pk+ 1/2  —  400  hPa,  followed  by  a  fairly  sharp  transi¬ 
tion  to  the  LG-coordinate  limit  of  r  — >  2  as  y  — >  yk  +1/2, 


4.  New  hybrid  coordinates 

The  NOGHYB  coordinate  developed  for  NOGAPS- 
ALPHA  was  chosen  to  have  the  following: 

1)  a  lowermost  SAL  model  layer  (rtT  =  1); 

2)  cr-like  properties  ( dBldy  constant)  over  the  —1000- 
400-hPa  range  (0.4  £  fj  S  1). 

Since  the  analytical  form  in  (30)  originally  designed 
for  NOGHYB  proves  general  enough  to  reproduce  all 
the  coordinates  in  section  3  (see  Table  1),  we  use  it  here 
to  develop  two  new  hybrid  coordinates  that  remove 
NOGHYB  restrictions  (1)  and  (2)  above. 
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Fig.  7.  NOGHYB  coordinate  response  of  layer  pressure  thick¬ 
nesses  A pk  vs  pk  (solid)  relative  to  the  L84  profile  A pk  at  ps  =  p0 
(dotted).  Gray  line  makes  pkp+1,2  86.7  hPa  (kp  =  58). 

a.  Relaxing  the  lower  boundary  limit:  Coordinate 

NEWHYB1 

Our  first  new  coordinate  (NEWHYB1)  elimi¬ 
nates  criterion  (1)  while  retaining  criterion  (2).  We 
achieve  this  by  setting  rp  =  2.2,  ra  =  1.2,  S  =  5,  and 
k a  =  0,  the  latter  to  avoid  thickness  discontinuities 
across  the  L  —  +  1/2  interface.  The  dashed  curve  in 

Fig.  8  plots  dBldj]  for  this  coordinate.  It  reveals  that  the 
limiting  surface  index  of  rL+la  =  ra.=  1.2  yields  a  sur¬ 
face  layer  that  is  now  neither  a  cr  nor  a  SAL  surface, 
so  that  criterion  (1)  no  longer  holds.  However,  dBldr\ 
remains  roughly  constant  from  0.4  s  fj  s  1,  satisfying 
criterion  (2).  The  rp  =  2.2  limit  provides  a  slightly  smooth¬ 
er  transition  to  isobaric  layers  than  the  original  rp  =  2, 
since  now  d2B!drf  — >  0  as  t)  — >  rj/,  +1/2  from  be¬ 
low.  Recall  from  (21)  that  vanishing  d2B/d fj2  im¬ 
plies  local  shape  preservation  of  the  thickness  profiles. 


Shape  factors  s(fjk,ps)  in  Fig.  9  retain  near-perfect 
shape  preservation  from  —1000  to  400  hPa  via  a  quasi¬ 
constant  s(r]k,ps)  profile,  followed  by  a  smooth  transi¬ 
tion  to  isobaric  layers  thereafter.  The  corresponding 
pressure  thickness  responses  A pk  in  Fig.  10  remain 
similar  in  shape  to  the  original  profile  over  high  terrain. 
The  shape  preservation  properties  of  this  new  coordi¬ 
nate  are  evidently  superior  to  those  of  NOGHYB  (cf. 
Fig.  7). 

b.  Preservation  of  thickness  profile  shapes  at  all 
model  levels:  Coordinate  NEWHYB2 

Here  we  also  relax  criterion  (2)  and  adjust  param¬ 
eters  to  attain  a  hybrid  coordinate  with  the  best  shape 
preservation  properties  not  just  over  the  —1000-400- 
hPa  range,  but  on  all  the  terrain-following  model  levels 
%+ 1/2  <*»<!• 

This  NEWHYB2  coordinate  is  identical  to 
NEWHYB1,  except  we  now  set  r„  =  1.35  (see  Table  1). 
This  coordinate’s  dBldr\  profile  is  plotted  as  the  dot- 
dashed  curve  in  Fig.  8.  The  resulting  shape  profiles 
s(v k>Ps)  in  Fig-  11  now  increase  gradually  over  the 
— 1000-400-hPa  range  and  thus  do  not  rigidly  preserve 
thickness  profile  shapes.  However,  when  applied  over 
the  whole  profile,  the  net  effects  of  the  new  coordinate 
in  preserving  the  smooth  nature  of  the  thickness  pro¬ 
files  are  superior.  Over  the  Himalayas,  for  example,  the 
NEWHYB2  coordinate  yields  a  compressed  thickness 
profile  A pk  in  Fig.  12d  that  better  preserves  the  original 
thickness  profile  shape  than  the  NEWHYB1  response 
in  Fig.  lOd.  We  have  tested  the  two  new  coordinates 
with  a  range  of  other  reference  thickness  profiles  like 
those  shown  in  Fig.  la,  and  have  found  in  each  case 
that,  of  the  three  NOGAPS-ALPHA  hybrids  in  Table 
1,  the  NEWHYB2  coordinate  consistently  exhibits  the 
best  overall  shape  preservation  over  high  terrain. 

Even  better  shape  preservation  is  obtained  by  in¬ 
creasing  r„  still  further:  for  example,  r,r  =  1.5  yields 
excellent  shape  preservation  properties.  However, 
these  larger  index  values  in  turn  yield  larger  values  of 
(Ps) min- 1 h  li s  making  the  hybrid  coordinate  more  prone 
to  violation  of  monotonicity  over  high  terrain.  Since 
NOGAPS-ALPHA  is  an  NWP  prototype  typically  run 


Table  1. 

The  different  vertical  coordinates  as  implemented  for  NOGAPS-ALPHA  L84.  For  all  these  coordinates  k a 
entries  for  S  indicate  that  its  value  is  not  relevant  in  setting  the  coordinate’s  behavior  since  rp  =  r,r 
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dB/cfx\ 

Fig.  8.  Profiles  of  dBldi )  for  the  new  hybrid  coordinates  dis¬ 
cussed  in  section  4,  with  the  NOGFIYB,  SAL,  and  cr  values  shown 
for  reference  (see  Table  1). 

at  high  horizontal  resolution  (see  Eckermann  et  al. 
2006),  we  settled  on  ra  =  1.35  as  a  reasonable  trade-off 
between  shape  preservation  and  robust  monotonicity. 


20  40 

A Pk  (hPa) 


60 


20  40 

Aft  (hPa) 


60 


(d)  90°E:  ps  =  469hPa 
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Fig.  10.  NEWHYB1  coordinate  response  of  layer  pressure 
thicknesses  A pk  vs  pk  (solid)  relative  to  the  L84  profile  A pk  at 
Ps  =  Po  (dotted). 


5.  Errors  in  pressure  gradient  force  over  steep 
orography 

In  cr-coordinate  models  the  pressure  gradient  force 
(PGF)  is  the  sum  of  two  terms.  Over  steeply  sloping 


terrain,  large  hydrostatic  corrections  apply  to  both 
terms,  making  their  values  very  nearly  equal  and  oppo¬ 
site.  The  small  PGF  residual  therefore  becomes  suscep¬ 
tible  to  the  combined  discretization  and  interpolation 
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Fig.  9.  NEWF1YB1  coordinate  shape  factors  s(rfk ,  ps )  for  L84 
hybrid  levels  in  Fig.  2f  at  the  indicated  longitudes  with  progres¬ 
sively  higher  terrain.  Thick  gray  curve  marks  pk  +1/2  86.7  hPa 

dP  58). 


Fig.  11.  NEWF1YB2  coordinate  shape  factors  s(fik,  ps)  for  L84 
hybrid  levels  in  Fig.  2g  at  the  indicated  longitudes  with  progres¬ 
sively  higher  terrain.  Thick  gray  curve  marks  pk  +1/2  -=  86.7  hPa 
dp  =  58). 
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Fig.  12.  NEWHYB2  coordinate  response  of  layer  pressure 
thicknesses  A pk  vs  pk  (solid)  relative  to  the  L84  profile  A pk  at 
Ps  =  Po  (dotted). 


errors  in  each  of  these  larger  magnitude  terms,  as  is  well 
known  (e.g.,  Corby  et  al.  1972;  Gary  1973;  Phillips  1974; 
Arakawa  and  Suarez  1983). 

Simmons  and  Burridge  (1981)  found  significant  re¬ 
ductions  in  stratospheric  PGF  errors  using  hybrid  co¬ 
ordinates.  To  facilitate  comparisons  with  their  results, 
here  we  evaluate  their  discretized  error  function  Ek{ps) 
using  their  vertical  temperature  profile  (Fig.  13a;  see 
also  Fig.  2  of  Simmons  and  Strafing  1981)  and  the  same 
choice  of  a  10%  surface  pressure  change  over  a  100-km 
horizontal  distance.  Figures  13b, c  show  resulting 
Ek(ps )  profiles  for  the  hybrid  coordinates  in  Table  1 
using  our  L84  model  layers.  As  in  Simmons  and  Bur¬ 
ridge  (1981),  these  values  are  plotted  as  geostrophic 
wind  speed  errors  using  a  nominal  Coriolis  parameter 
/  =  10-4  s_1.  We  show  results  for  ps  =  800  hPa.  The 
abrupt  increase  in  stratospheric  temperatures  leads  to 
large  stratospheric  PGF  errors  using  the  <j  coordinate. 
All  the  hybrid  coordinates  reduce  these  stratospheric 
errors  while  maintaining  fairly  small  tropospheric 
errors.  Figure  13c  shows  that  the  NEWHYB1  and 
NEWHYB2  coordinates  yield  the  smallest  strato¬ 
spheric  PGF  errors  of  all  the  hybrid  coordinates  tested. 

To  investigate  the  generality  of  these  findings,  we 
allowed  the  analytical  temperature  profile  to  vary 
somewhat  realistically.  First  we  added  normally  distrib¬ 
uted  random  offsets  to  the  sea  level  and  tropopause 


temperatures  of  3-K  standard  deviation,  to  stratopause 
temperatures  of  5-K  standard  deviation,  and  random 
vertical  offsets  to  the  tropopause  pressure  height  of 
1-km  standard  deviation.  Then,  to  each  temperature 
profile  j  we  superimposed  a  sinusoidal  perturbation  of 
the  following  form: 


Tj(Z)  =  Tpeak(Z)  sin 


w*(z)  dz  +  <pj 


Tpeak(Z)  =  IK  exp  I  - 


/  Z  —  10  km 


2HP 


•M-Z)  - 


2tt 

mJZ) 


=  2  km 


expl - 


/  Z  —  10  km 


2  Hp 


(32) 

(33) 

(34) 


These  equations  mimic  the  vertical  variation  of  the 
energy-containing  gravity  wave  scales  as  defined  by  ca¬ 
nonical  vertical  wavenumber  spectra  (see  Figs.  2-3  of 
Smith  et  al.  1987),  with  parameters  settings  in  (33)  and 
(34)  based  on  spectral  analyses  of  long-term  high- 
resolution  radiosonde  temperatures  (e.g.,  Allen  and 
Vincent  1995).  We  set  HE  =  2.3//  16  km  following 

Fritts  and  Alexander  (2003),  and  assign  phases  ran¬ 
domly  between  0°  and  360°.  Ten  perturbed  tempera¬ 
ture  profiles  resulting  from  this  algorithm  are  plotted  in 
Fig.  13d  to  convey  the  scale  and  range  of  the  imposed 
variability. 

The  root-mean  square  errors  (RMSEs)  calculated 
from  500  such  profiles  using  ps  =  800  hPa  are  plotted  in 
Figs.  13e,f  for  the  various  hybrid  coordinates.  The  over¬ 
all  trends  seen  in  the  top  panels  for  the  original  tem¬ 
perature  profile  persist.  In  particular,  the  NEWHYB1 
and  NEWHYB2  coordinates  yield  the  smallest  strato¬ 
spheric  PGF  errors. 


6.  Tests  in  a  global  model 

In  addition  to  PGF  errors,  hybrid  coordinates  can 
alleviate  other  sources  of  cr-coordinate  errors  that  are 
not  easily  gauged  by  offline  tests,  such  as  errors  in  simu¬ 
lating  vertical  advection  or  sharp  horizontal  gradients 
(Konor  and  Arakawa  1997).  Thus,  here  we  implement 
a  subset  of  the  hybrid  vertical  coordinates  introduced  in 
sections  3  and  4  into  a  global  spectral  model  and  test 
their  performance  using  direct  model  integrations.  We 
use  the  forecast  model  component  of  NOGAPS- 
ALPHA,  the  high-altitude  NOGAPS  prototype,  which 
is  hydrostatic,  Eulerian,  spectral  in  the  horizontal  and 
finite  difference  in  the  vertical  based  on  the  L-grid  dis¬ 
cretization  of  Arakawa  and  Suarez  (1983)  (see  also 
Hogan  and  Rosmond  1991).  A  vertical  coordinate  of 
Simmons  and  Striifing’s  (1983)  form  in  (6)  is  built  into 
the  dynamical  core. 
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Fig.  13.  For  the  idealized  temperature  profile  in  (a),  (b)-(c)  plot  L84  geostrophic  wind  errors  for  the  indicated 
coordinates  based  on  the  specific  vertical  discretization  of  Simmons  and  Burridge  (1981):  see  their  Eq.  (6.1).  (e), 
(f)  Geostrophic  wind  RMSEs  from  500  similar  calculations  applied  to  temperature  profiles  with  stochastic  vari¬ 
ability  added,  as  described  in  section  5:  (d)  10  representative  examples  of  the  temperature  profiles.  All  calculations 
used  a  surface  pressure  ps  =  800  hPa. 


a.  Idealized  simulations 

The  adiabatic  simulations  using  analytical  initial  con¬ 
ditions  that  are  routinely  used  to  test  global  models 
(e.g.,  Giraldo  and  Rosmond  2004)  typically  cither  omit 
terrain  (Staniforth  and  White  2007)  or  prescribe  spe¬ 
cific  planetary-scale  terrain  variations  (Jablonowski 
and  Williamson  2006).  To  test  hybrid  coordinates,  we 
require  states  that  support  more  realistic  variable  ter¬ 
rain.  Static  atmospheres  in  which  temperature  only 
varies  vertically  permit  arbitrary  terrain  since,  given 
an  analytical  temperature  profile,  hydrostatic  integra¬ 
tion  yields  corresponding  analytical  expressions  for  the 
required  surface  temperature  and  pressure  as  a  func¬ 
tion  of  terrain  height  needed  to  initialize  this  state. 
Static  states  also  allow  dynamical  errors  to  be  diag¬ 
nosed  directly  from  the  simulated  wind,  divergence, 
and  vorticity  fields.  T79L84  NOGAPS-ALPHA  ex¬ 
periments  using  the  linearly  decreasing  temperature 
profile  in  Giraldo  and  Rosmond  (2004)  maintained  a 


static  state  over  15  days  using  both  ideal  and  realistic 
terrain.  We  then  used  the  COESA  (1976)  temperature 
profile,  which  includes  realistic  middle-atmosphere 
variations  as  in  Fig.  13a.  We  show  results  using  the 
latter  temperature  initialization  and  realistic  T79  orog¬ 
raphy  only,  as  they  are  very  similar  to  those  using  ide¬ 
alized  terrain. 

Figure  14  plots  global  RMSEs  of  relative  vorticity 
and  divergence  versus  forecast  hour  at  500,  100,  and  10 
hPa.  These  errors,  while  still  small,  are  about  two  to 
three  orders  of  magnitude  larger  than  those  using  the 
linear  temperature  profile.  The  vorticity  errors  mani¬ 
fest  over  the  terrain  as  vertically  deep  structures  of 
short  horizontal  scale,  consistent  with  the  PGF  errors 
discussed  in  section  5.  Like  the  offline  PGF  estimates 
in  Fig.  13,  the  vorticity  errors  are  largest  for  the  cr  co¬ 
ordinate  at  all  levels,  and  are  one  to  two  orders  of 
magnitude  smaller  when  using  the  NEWFIYB1  or 
NEWHYB2  coordinate. 

Divergence  errors  in  Fig.  14  are  also  largest  using  cr 
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Fig.  14.  Variation  with  forecast  time  of  global  RMSEs  in  (top)  relative  vorticity  and  (bottom)  divergence  at  (from  left  to  right)  500, 
100,  and  10  hPa  for  idealized  T79L84  NOGAPS-ALPHA  experiments  using  different  vertical  coordinates:  a  (solid),  SAL  (dotted), 
NOGHYB  (dashed),  NEWHYB1  (dot-dashed),  and  NEWHYB2  (solid-dashed). 


coordinates.  They  aggregate  over  high  steep  terrain 
such  as  the  Andes  and  Himalayas  as  short  shallow 
structures,  very  similar  to  the  spurious  stratospheric  di¬ 
vergence  structures  reported  by  Trenberth  and  Stepa- 
niak  (2002)  over  high-terrain  in  cr-based  reanalysis 
fields.  The  results  in  Fig.  14  support  the  contention  of 
Trenberth  and  Stepaniak  (2002)  that  these  anomalies 
should  be  substantially  reduced  by  using  hybrid  coor¬ 
dinates.  Again,  the  error  suppression  is  greatest  at  all 
altitudes  when  using  the  NEWHYB1  or  NEWHYB2 
coordinate. 

b.  Stability  of  the  semi-implicit  time  scheme 

NOGAPS  achieves  longer  model  time  steps  by  treat¬ 
ing  gravity  wave  propagation  semi-implicitly.  The  algo¬ 
rithm,  described  in  section  7.2  of  Hogan  et  al.  (1991), 
closely  follows  the  discretized  three-time-level  (leap¬ 
frog)  implementation  of  the  Robert  et  al.  (1972) 
method  by  Hoskins  and  Simmons  (1975),  as  subse¬ 
quently  generalized  to  hybrid  vertical  coordinates  by 
Simmons  and  Striding  (1981).  Following  Hoskins  and 
Simmons  (1975),  the  algorithm  linearizes  around  a  ref¬ 
erence  temperature  profile  Tr{ rj).  Based  on  the  stability 
analysis  of  Simmons  et  al.  (1978),  NOGAPS  uses  an 
isothermal  reference  temperature  profile  of  T,.( -rj)  = 
300  K  (Hogan  et  al.  1991),  the  same  choice  used  in  the 
hybrid  coordinate  version  of  the  ECMWF  model  (Sim¬ 
mons  et  al.  1989). 


Extension  of  this  scheme  to  hybrid  coordinates  in¬ 
volves  linearization  about  a  reference  surface  pressure, 
pr,  as  well  (Simmons  and  Burridge  1981).  Since 
NOGAPS  uses  cr  levels  operationally,  the  choice  for  pr 
has  not  been  critical  to  the  stability  of  the  scheme  to 
date  (see,  e.g..  Table  3  of  Simmons  and  Burridge  1981). 
Hogan  et  al.  (1991)  adopted  pr  =  600  hPa,  based  on 
calculations  suggesting  a  small  stability  advantage  in 
using  pr  values  near  the  minimum  anticipated  terrain 
pressures  in  the  model.  In  testing  various  types  of  hy¬ 
brid  coordinates,  however,  Simmons  and  Strafing 
(1981)  found  that  the  stability  of  the  semi-implicit 
scheme  became  sensitive  to  the  value  of  pr. 

Simmons  et  al.  (1978)  and  Simmons  and  Strafing 
(1981)  modeled  the  stability  of  their  semi-implicit  time 
scheme  offline  by  using  the  idealized  vertical  tempera¬ 
ture  profile  in  Fig.  13a,  thereby  simplifying  the  problem 
to  a  matrix  eigenvalue  equation  that  they  solved  nu¬ 
merically  (see  appendix  A  of  Simmons  et  al.  1989). 
Here  we  test  the  stability  of  the  NOGAPS  semi-implicit 
time  scheme  by  running  NOGAPS-ALPHA  and  moni¬ 
toring  its  output.  We  favor  this  “brute  force”  approach 
over  simpler  offline  eigenvalue  computations  since  the 
latter  provides  only  broad  guidance,  tends  to  overesti¬ 
mate  the  scheme’s  stability  relative  to  full  model  imple¬ 
mentations  (Simmons  et  al.  1978),  and  sometimes  yields 
results  counter  to  the  experience  from  model  runs 
(Simmons  et  al.  1989). 
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Table  2.  Standard  NOGAPS  model  time  steps  At  (s)  and  their  corresponding  integer  fractions  of  1  h. 


At  (s) 

1800 

1200 

900 

720 

600 

450 

400 

360 

300 

240 

200 

180 

150 

120 

100 

1  hi  At 

2 

3 

4 

5 

6 

8 

9 

10 

12 

15 

18 

20 

24 

30 

36 

We  test  model  stability  by  running  the  same  T79L84 
experiment  in  every  case.  The  run  is  initialized  with 
archived  analysis  fields  at  0000  UTC  14  January  2003 
and  run  out  for  46  days  to  0000  UTC  1  March.  In  com¬ 
mon  with  other  models,  NOGAPS-ALPHA  generally 
updates  its  radiative  heating  and  cooling  rates  only  ev¬ 
ery  1  or  2  h.  In  order  that  radiation  and  model  time 
steps  coincide,  only  model  time  steps  At  for  which  1  hi  At 
is  an  integer  are  permitted.  Table  2  lists  the  allowed 
values,  and  we  run  with  those  successively  smaller  time 
steps  until  a  stable  46-day  model  integration  is 
achieved.  A  stable  run  is  identified  if  all  model  fields 
remain  finite  over  the  entire  46  days.  The  resulting 
stable  At  values  are  listed  in  Table  3  as  a  function  of 
hybrid  coordinate  and  pr.  To  assess  overall  perfor¬ 
mance  objectively,  Table  3  also  compiles  scores  by  add¬ 
ing  At  values  in  each  row  and  column. 

Table  3  reveals  instability  for  every  coordinate  when 
pr  =  500  hPa.  Similar  uniform  instability  at  pr  =  500 
hPa  was  noted  in  the  eigenvalue  solutions  of  Simmons 
and  Strafing  (1981)  and  Simmons  and  Burridge  (1981). 
At  600  hPa  (the  current  default  NOGAPS  setting),  the 
cr,  SAL,  CMAM,  and  NOGHYB  coordinates  yield 
stable  integrations  at  At  ~  600-720  s,  whereas  the  LG, 
NEWHYB1,  and  NEWHYB2  coordinates  only  yield 
stable  integrations  at  much  smaller  time  steps.  Thus, 
the  transition  from  cr  to  progressively  more  aggressive 
hybrid  coordinates  gradually  reduces  the  low-end  range 
of  pr  values  for  which  the  semi-implicit  time  scheme  is 
stable  at  acceptably  long  model  time  steps.  This  is  re¬ 
flected  in  the  column  scores  along  the  bottom  row  of 


Table  3,  which  are  large  for  u  and  small  for  NEWHYB2 
and  LG.  Interestingly,  the  highest  column  score  occurs 
for  the  NOGHYB  coordinate.  Since  NOGHYB  was 
designed  as  a  first  hybrid  replacement  for  cr  levels  in 
NOGAPS-ALPHA  (see  section  3d),  the  results  in 
Table  1  indicate  that  it  offers  no  stability  penalty  (in¬ 
deed,  slight  advantages)  relative  to  cr  at  any  pr  >  600 
hPa.  This  objectively  confirms  anecdotal  experience 
with  this  coordinate  in  NOGAPS-ALPHA,  which, 
since  it  was  implemented,  has  yielded  stable  integra¬ 
tions  in  a  variety  of  short-  and  long-term  simulations  at 
the  default  NOGAPS  setting  of  pr  =  600  hPa  (Allen  et 
al.  2006;  McCormack  et  al.  2006;  Eckermann  et  al. 
2007a). 

Furthermore,  for  pr  >  700  hPa  the  NEWHYB1  and 
NEWHYB2  coordinates  also  exhibit  equal  or  superior 
stability  to  cr.  The  highest  row  scores  in  the  far-right 
column  of  Table  3  occur  for  pr  =  800  hPa,  indicating 
that  this  value  yields  the  highest  overall  stability  of  the 
semi-implicit  time  scheme  for  all  the  hybrid  coordinates 
tested.  This  is  a  robust  result,  evident  in  Table  3  for  cr 
levels  all  the  way  through  to  the  most  aggressive  hybrid 
coordinates  like  LG  and  NEWHYB2.  These  findings 
are  consistent  with  the  choice  of  pr  =  800  hPa  in  the 
corresponding  ECMWF  scheme  (Simmons  et  al.  1989). 

The  outlier  is  the  SAL  coordinate,  which  is  the  only 
coordinate  tested  that  does  not  exhibit  a  regular  trend 
in  stability  as  a  function  of  pr  in  Table  3.  As  discussed 
in  section  3b,  the  large  abrupt  discontinuities  in  layer 
thicknesses  at  pk  +1/2  (see  Figs.  4  and  5)  lead  to  inac¬ 
curate  vertical  finite-difference  calculations  that  reduce 


Table  3.  Longest  model  time  steps  A t  (s;  see  Table  2)  yielding  stable  integration  of  the  46-day  NOGAPS-ALPHA  T79L84  hindcast, 
as  a  function  of  vertical  coordinate  (see  Table  1 )  and  pr.  Total  “scores”  are  given  in  the  bottom  row  and  right  column  by  adding  At  values 
along  the  indicated  column  and  row,  respectively.  Table  entry  “U”  denotes  unstable  integrations  for  all  model  time  steps  listed  in  Table 
2  and  is  assigned  a  score  of  zero. 


Vertical  coordinates 


Pr 

( 7 

SAL 

CMAM 

NOGHYB 

NEWHYB1 

NEWHYB2 

LG 

Tot  score 

500  hPa 

u 

U 

U 

U 

U 

U 

U 

0 

550  hPa 

600 

240 

100 

100 

U 

u 

U 

1040 

600  hPa 

600 

720 

600 

720 

200 

150 

100 

3090 

700  hPa 

600 

720 

600 

720 

720 

600 

600 

4560 

800  hPa 

720 

600 

600 

720 

720 

720 

720 

4800 

900  hPa 

600 

450 

600 

720 
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600 

600 

4290 

1000  hPa 

600 

720 

600 

720 

720 

600 

600 

4560 

1100  hPa 

600 

600 

600 

720 

600 

600 

600 

4320 

Tot  score 

4320 

4050 

3700 

4420 
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Fig.  15.  Northern  Hemisphere  50-hPa  geopotential  heights  at  0000  UTC  (top)  13,  (middle)  15,  and  (bottom)  19  Jan  2003.  (left) 
NOGAPS  reanalysis  fields.  Remaining  plots  show  NOGAPS-ALPHA  T239L84  +  2-,  +4-,  and  +8-day  forecasts  initialized  with 
NOGAPS/GMAO  reanalysis  fields  at  0000  UTC  11  Jan  2003  using  NOGHYB,  NEWHYB1,  NEWHYB2,  and  a  vertical  coordinates 
(columns  2-5,  respectively).  Contour  interval  is  160  m  with  the  range  19  440-20  080  m  shaded  gray. 


the  stability  of  the  semi-implicit  time  scheme  (see  sec¬ 
tion  7  of  Simmons  and  Burridge  1981).  Indeed,  Sim¬ 
mons  and  Striifing  (1983)  found  that  the  ECMWF 
model  was  more  prone  to  instability  with  the  SAL  co¬ 
ordinate  than  with  any  other  coordinate  they  tested, 
consistent  with  the  results  in  Table  3. 

Thus,  changing  pr  from  600  to  800  hPa  in  NOGAPS- 
ALPHA  provides  optimal  model  stability  over  the 
broad  range  of  hybrid  coordinates  tested  in  Table  3. 
More  specifically,  this  change  allows  the  NEWHYB1 
and  NEWHYB2  coordinates  to  be  implemented  and 
tested  in  NOGAPS-ALPHA  without  incurring  any  sta¬ 
bility  penalty  relative  to  cr. 

c.  Forecast  skill 

To  assess  the  overall  impact  of  each  hybrid  coordi¬ 
nate  in  NOGAPS-ALPHA,  we  performed  T239L84 
forecasts  using  each  coordinate,  initialized  at  0000  UTC 
for  each  day  of  the  period  6-30  January  2003.  We  chose 
this  period  because  of  the  range  of  stratospheric  dy¬ 
namics  that  occurred.  The  early  parts  of  January  2003 
contained  strong  stratospheric  vortex  winds  (see  Fig.  8 
of  Eckermann  et  al.  2007b),  which  present  well-known 


challenges  for  model  stability  (e.g.,  Simmons  et  al. 
1989).  In  mid-January  there  was  a  major  stratospheric 
warming  (McCormack  et  al.  2004;  Blum  et  al.  2006)  of 
the  split-vortex  type  (Charlton  and  Polvani  2007),  like 
the  February  1979  warming  focused  upon  in  the  hybrid 
coordinate  tests  of  the  ECMWF  model  by  Simmons 
and  Striifing  (1983)  and  the  February  1994  warming  in 
the  forecasting  study  of  Lahoz  (1999).  Since  warmings 
present  the  greatest  stratospheric  forecasting  and  sta¬ 
bility  challenges  for  NWP  models  (Waugh  et  al.  1998; 
Lahoz  1999;  Simmons  et  al.  2005),  this  period  provides 
a  stringent  test  of  model  stability  and  forecast  skill. 

As  both  initial  states  and  verifying  analysis  for  the 
forecasts,  we  use  NOGAPS  reanalyses  (1000-10  hPa) 
from  the  T239L30  operational  (cr  coordinate)  sys¬ 
tem  using  the  Naval  Research  Laboratory  (NRL)  At¬ 
mospheric  Variational  Data  Assimilation  System 
(NAVDAS),  the  National  Aeronautics  and  Space  Ad¬ 
ministration’s  (NASA’s)  Global  Modeling  and  Assim¬ 
ilation  Office  (GMAO)  analyses  (10-0.4  hPa),  and  a 
relaxation  to  climatology  above  the  0.4-hPa  level.  De¬ 
tails  on  the  blending  and  model  initialization  of  these 
fields  are  given  in  section  3.1.2  of  Eckermann  et  al. 
(2006). 
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Fig.  16.  (top  two  rows)  Mean  geopotential  height  ACCs  and  (bottom  row)  horizontal  wind 
RMSEs  vs  forecast  hour  (6-h  resolution)  for  T239L84  NOGAPS-ALPHA  forecasts  from  6  to 
30  Jan  2003.  (top  insets)  Results  for  different  vertical  coordinates:  cr  (solid),  SAL  (dotted), 
NOGHYB  (solid-dashed),  NEWHYB1  (dot-dashed),  and  NEWHYB2  (dashed),  (bottom 
insets)  Differences  between  the  latter  four  curves  and  the  cr  curve,  with  signs  that  signify  better 
or  worse  skill  than  the  cr  forecast  highlighted  in  each  case  in  gray.  Results  are  shown  for  (a) 
500  hPa  20°-80°S,  (b)  500  hPa  20°-80°N,  (c)  100  hPa  20°-80°S,  (d)  100  hPa  20°-80°N,  (e)  50 
hPa  20°-80°N,  and  (f)  5  hPa  20°-80°N. 


Each  forecast  ran  out  10  days,  with  output  saved  ev¬ 
ery  6  h.  The  forecasts  all  used  p,.  =  800  hPa  based  on  the 
results  in  Table  3.  As  an  example,  Fig.  15  shows  50-hPa 
geopotential  heights  from  +2-,  +4-,  and  -l-8-day  fore¬ 
casts  initialized  on  11  January  using  different  vertical 
coordinates,  with  the  NOGAPS/GMAO  reanalysis 
fields  shown  in  the  left  column.  The  +  8-day  forecasts 
capture  the  onset  of  a  split-vortex  warming.  To  assess 
the  skill  of  all  the  forecasts  over  this  25-day  period 
objectively,  we  computed  mean  errors  (bias),  mean  ab¬ 


solute  errors,  RMSEs,  and  anomaly  correlation  coeffi¬ 
cients  (ACCs)  using  standard  area-weighted  formulas 
(see  section  2  of  Waugh  et  al.  1998). 

Figure  16  plots  time-mean  geopotential  height  ACCs 
for  these  forecasts  at  500  (top  row)  and  100  hPa  (middle 
row)  in  the  20°-80°  latitude  band  in  each  hemisphere. 
The  geopotential  height  climatology  used  in  these  cal¬ 
culations  was  the  15-yr  January  mean  of  the  NCEP 
Department  of  Energy  (DoE)  Second  Atmospheric 
Model  Intercomparison  Project  (AMIP  II)  reanalysis 
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(Kanamitsu  et  al.  2002).  There  is  very  little  differ¬ 
ence  among  the  ACCs  of  the  various  coordinate  fore¬ 
casts,  even  beyond  5  days.  Similarly  small  (neutral) 
model  impacts  of  hybrid  coordinate  forecasts  relative  to 
cr  forecasts  were  noted  by  Simmons  and  Strafing 
(1983).  The  small  differences  are  expanded  in  the  bot¬ 
tom  panel  of  each  plot  by  plotting  the  difference  in  each 
hybrid  coordinate’s  ACC  from  the  ACC  of  the  cr- 
coordinate  forecasts.  These  plots  reveal  small  positive 
impacts  in  the  Southern  Hemisphere  and  small  negative 
impacts  in  the  Northern  Hemisphere.  Although  all 
these  trends  are  very  weak,  close  inspection  shows  that 
the  NEWHYB1  and  NEWHYB2  coordinate  forecasts 
have  small  systematic  positive  impacts  relative  to  the 
other  hybrid  coordinate  forecasts. 

Figures  16e,f  show  stratospheric  horizontal  wind 
RMSEs  in  the  20°-80°N  band.  Since  the  error  analysis 
in  sections  5  and  6a  indicates  that  differences  should  be 
largest  over  high  variable  terrain,  these  RMSE  calcula¬ 
tions  used  winds  only  from  those  model  grid  points 
where  the  underlying  terrain  was  >1000  m.  While 
RMSE  differences  are  again  very  small,  overall  the  hy¬ 
brid  coordinate  RMSEs  are  systematically  lower  than 
the  cr  reference.  This  appears  to  be  consistent  with  the 
idealized  PGF  error  estimates  in  Figs.  13  and  14. 

7.  Summary  and  discussion 

Perhaps  the  central  result  of  this  study  is  clarification 
of  the  role  of  dBldr)  in  defining  salient  characteristics  of 
hybrid  coordinates  of  the  Simmons  and  Strafing  (1983) 
form  in  (6).  In  diagnosing  coordinate  properties  using 
(7/l/c/f)-based  tests  of  monotonicity  and  preservation  of 
the  smoothness  and  shape  of  layer  thickness  profiles,  it 
became  clear  that  there  are  no  obvious  optimal  choices 
for  this  profile:  each  choice  has  specific  advantages  and 
disadvantages. 

For  example,  robust  shape  preservation  and  mono¬ 
tonicity  competitive  with  cr  levels  suggests  a  SAL  or 
NOGHYB  coordinate  with  constant  or  near-constant 
tropospheric  dBldr]  values.  Yet  this  property  then  re¬ 
quires  a  very  sharp  transition  in  the  lower  stratosphere 
to  yield  isobaric  model  layers  in  the  middle  atmosphere 
(Fig.  3),  which  distorts  layer  thickness  profiles  over  high 
terrain  in  the  upper  troposphere  and  lower  stratosphere 
(see  Figs.  4,  5,  and  7).  This,  then,  would  seem  to  argue 
for  hybrid  coordinates  incorporating  a  smoother  varia¬ 
tion  in  dBldr]  with  height,  such  as  the  LG  coordinate. 
Yet  for  these  coordinates,  the  integral  constraint  on 
dBldr)  from  the  boundary  conditions  (10), 

1  dB 

—  (!]')  dr]' =B(  1)  =  1,  (35) 

Vkp+V2  d  V 


implies  necessarily  larger  maximum  dBldr]  values  that 
render  these  coordinates  nonmonotonic  at  larger 
( /?s)min  values.  The  final  choice  involves  a  trade-off  be¬ 
tween  these  and  other  competing  concerns,  and  thus 
will  vary  from  model  to  model.  For  instance,  low- 
resolution  climate  models  need  worry  less  about  the 
monotonicity  criterion  implied  by  maximum  dBldr]  val¬ 
ues  than  high-resolution  NWP  systems. 

In  choosing  a  coordinate  it  is  instructive  to  review  the 
choices  made  in  other  models.  Some  analytical  choices 
were  reviewed  in  section  3.  More  commonly,  however, 
modeling  studies  simply  tabulate  their  discretized 
Ak+ V2  and  Bk+ln  values.  Direct  inspection  of  these  val¬ 
ues  is  not  particularly  illuminating  and,  thus,  is  of  lim¬ 
ited  guidance  to  other  modelers.  However,  dBldr]  can 
readily  be  evaluated  from  such  tabulations,  which,  we 
have  argued,  provides  all  the  salient  information  about 
the  nature  of  the  hybrid  coordinate.  With  this  under¬ 
standing,  we  can  compute  and  plot  dBldr]  values  for 
different  models. 

Figure  17  plots  dBldr]  profiles  computed  from  quoted 
B k +1/2  coefficients  used  in  some  well-known  global  cli¬ 
mate  models,  with  the  source  of  these  values  cited  in 
each  example  in  the  figure  caption.  Unsurprisingly, 
based  on  our  previous  discussion,  hybrid  coordinate 
choices  vary  widely  among  the  models.  For  example. 
Fig.  17a  immediately  reveals  that  the  Community  At¬ 
mosphere  Model  (CAM)  and  Whole  Atmosphere 
Community  Climate  Model  (WACCM)  use  a  SAL  co¬ 
ordinate  (cf.  Fig.  3).  Additionally,  we  see  that  their  low¬ 
est  model  layer  is  set  to  pure  cr  ( dBldr]  =  1),  which 
produces  a  small  low-level  thickness  discontinuity,  in 
addition  to  the  standard  SAL  thickness  discontinuity  at 
Pkp+ 1/2- 

Similar  Bk+1/ 2  tabulations  for  NWP  systems  are 
harder  to  find  in  the  scientific  literature.  Figure  18 
plots  dBldr]  profiles  from  six  different  versions  of  the 
ECMWF  IFS,  obtained  from  Bk+m  tabulations  made 
publicly  available  on  the  ECMWF  Web  site,  and  those 
from  the  GMAO  Goddard  Earth  Observing  System 
(GEOS4)  and  the  German  Weather  Service  model  ob¬ 
tained  from  technical  reports.  The  ECMWF  results  re¬ 
veal  a  history  of  periodically  swapping  between  two 
different  types  of  hybrid  coordinate,  while  always  as- 
ymptoting  to  a  pure  cr  layer  at  the  surface.  The  German 
NWP  formulation  mirrors  one  of  these  two  ECMWF 
choices  (as  does  the  ECHAM5  formulation  in  Fig.  17e). 
The  pure  cr  layer  at  the  surface  of  the  GEOS4  model 
produces  a  sharp  discontinuity  in  layer  thicknesses  be¬ 
tween  the  two  lowest  layers  in  this  model. 

We  developed  and  analyzed  two  new  hybrid  coordi¬ 
nates  (NEWHYB1  and  NEWHYB2)  as  potential  re¬ 
placements  for  the  current  NOGHYB  coordinate  in 
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Fig.  17.  Profiles  of  dBidr\  used  in  selected  climate  models:  (a)  CAM3/WACCM  (Collins  et 
al.  2004);  (b)  GFDL  AM2  (Table  2  of  Anderson  et  al.  2004);  (c)  CCC  GCMII  (Table  1  of 
McFarlane  et  al.  1992);  (d)  CMAM  [Eq.  (1)  of  Polavarapu  et  al.  2005];  (e)  ECHAM5  (Table 
2.2  of  Roeckner  et  al.  2003);  and  (f)  ARPEGE  (Table  1  of  Deque  et  al.  1994). 


NOGAPS-ALPHA.  Offline  tests  suggested  these  new 
coordinates  offered  significant  improvements  relative 
to  NOGHYB  in  terms  of  shape  preservation  of  layer 
thickness  profiles  (cf.  Figs.  7, 10,  and  12)  and  reduction 
of  stratospheric  PGF  errors  (Figs.  13c, f).  While  ideal¬ 
ized  NOGAPS-ALPHA  experiments  confirmed  these 
error  reductions  (Fig.  14),  detailed  forecasting  experi¬ 
ments  indicated  essentially  neutral  (weak  positive)  im¬ 
pact  relative  to  NOGHYB  and  even  cr  forecasts.  Slight 
destabilization  of  the  semi-implicit  time  scheme  was 
rectified  by  changing  pr  from  600  to  800  hPa,  yielding 
stability  equal  to  or  better  than  cr.  Skill  scores  in  terms 
of  geopotential  height  ACCs  and  wind  RMSEs  all 
showed  marginal  differences,  with  small  reductions  in 
wind  errors  in  the  stratosphere.  Similar  neutral  or 
weakly  positive  findings  were  reported  by  Simmons  and 
Strafing  (1983),  suggesting  that  model  forecast  errors 
are  being  dominated  by  other  aspects  of  the  system, 
such  as  initial  conditions  and  physical  parameteriza- 
tions. 

Thus,  despite  offline  and  idealized  online  evidence  to 
the  contrary,  there  is  limited  practical  evidence  emerg¬ 


ing  from  short-term  forecast  integrations  to  forcefully 
recommend  one  hybrid  coordinate  choice  over  another. 
These  findings  in  some  sense  objectively  ratify  the  wide 
range  of  hybrid  coordinates  that  has  emerged  naturally 
among  the  current  generation  of  NWP  and  climate 
models  (Figs.  17  and  18).  The  one  clear  exception  to 
this  statement  is  the  SAL  coordinate,  which  has  ren¬ 
dered  both  NOGAPS-ALPHA  and  the  ECMWF 
model  of  Simmons  and  Strafing  (1983)  more  numeri¬ 
cally  unstable  than  any  other  vertical  coordinate,  and 
thus  seems  to  be  worth  avoiding.  The  following  quali¬ 
fying  remarks  also  pertain.  There  has  been  no  system¬ 
atic  investigation  of  the  impacts  of  hybrid  coordinate 
choice  on  either  extended-range  forecasts  or  climate 
integrations.  Furthermore,  the  short-term  forecasting 
tests  of  this  study  and  of  Simmons  and  Strafing  (1983) 
applied  only  to  the  forecast  model  component:  both 
studies  used  archived  analysis  fields  produced  by  a  cr- 
coordinate  model  for  forecast  initialization  and  verifi¬ 
cation.  The  inclusion  of  hybrid  coordinates  within  the 
full  forecast-assimilation  system  is  expected  to  yield 
further  improvements  in  skill  during  the  data  assimila- 


242 


MONTHLY  WEATHER  REVIEW 


Volume  137 


(a)  ECMWF  IFS  L 1 9  (b)  ECMWF  IFS  L3 1 


(g)  GMAO  GEOS4  L55 


(h)  German  NWP 


Fig.  18.  Profiles  of  dBldr]  used  in  selected  NWP  models.  Values  from  various  generations 
of  the  ECMWF  IFS  that  used  (a)  L19,  (b)  L31,  (c)  L40,  (d)  L50,  (e)  L60,  and  (f)  L91  layer 
thicknesses  (as  tabulated  online  at  http://www.ecmwf.int).  Results  for  (g)  GMAO  GEOS4  L55 
(Bloom  et  al.  2005)  and  (h)  German  Weather  Service  L31  Global-Modell  (GME)  (DWD 
2004). 


lion  phase  as  well  as  improving  the  overall  computa¬ 
tional  efficiency  of  the  system  (Simmons  and  Striifing 
1983;  Trenberth  and  Stepaniak  2002).  Given  the  large 
anomalies  in  stratospheric  cr-coordinate  analyses  re¬ 
ported  by  Trenberth  and  Stepaniak  (2002),  hybrid  co¬ 
ordinate  effects  on  data  assimilation  appear  to  be  worth 
studying  in  future  investigations. 

Based  on  the  consensus  of  the  tests  reported  here,  we 
have  recently  changed  the  default  hybrid  coordinate  in 
NOGAPS-ALPHA  from  NOGHYB  to  NEWHYB2, 
and  changed  the  default  pr  from  600  to  800  hPa. 


Acknowledgments.  This  research  was  supported  by 
the  Office  of  Naval  Research  and  the  Defense  Threat 
Reduction  Agency.  NCEP-DoE  reanalysis  2  data  were 
provided  by  the  NOAA/OAR/ESRL  PSD,  Boulder, 
CO  (see  online  at  http://www.cdc.noaa.gov).  Forecast 
runs  were  made  possible  by  computer  time  at  the  Air 
Force  Research  Laboratory  Major  Shared  Resource 
Center  provided  by  the  DoD  Fligh  Performance  Com¬ 
puting  Modernization  Program.  I  thank  Tim  Hogan  for 
generating  the  NOGAPS  reanalysis  fields  and  Frank 
Giraldo  for  providing  his  test  case  subroutines. 


January  2009 


ECKERMANN 


243 


APPENDIX 

Algorithm  for  Assigning  Ak+ ,/2  and  Bk+ V2 

The  starting  point  is  the  model’s  predefined  array  of 
layer  thicknesses  A pk  (see  Fig.  la).  Starting  at  the  top 
half-level  pV2  =  ptop,  these  thicknesses  are  integrated  to 
generate  interface  pressures: 

k 

Pk+ 1/2  =  Ptop  +  2  APf  (A1) 

i  I 

The  thicknesses  are  normalized  so  that  the  reference 


surface  pressure  ps  =  pL+1/2  in  (Al)  is  —1000  hPa.  The 
corresponding  half-level  rj  values  are  then 


Vk+l/2 


Pk+l/2  P  top 
P  L+l/2  Ptop 


(A2) 


The  algorithm  requires  specification  of  the  number 
of  isobaric  layers  at  the  top  of  the  model,  kp ,  and 
number  of  purely  cr-like  levels  at  the  bottom  of  the 
model,  ka¬ 
li  kp  +  0,  the  kp  isobaric  layers  are  set  by  assigning 
the  interface  coefficients  as 


Afc+l/2  —  Pk+l/2  —  P  top  + 
Bk+ 1/2  =  0 


/  Vk  +  l/2 
\  Vkp+l/2 


(Pk 


P  to 


1 . kp. 


(A3) 


Hybrid  cr-p  layers  are  set  by  the  following  relations  that  transition  from  a  pure  pressure  interface  at  kp  +  1/2 
to  a  pure  cr-likc  interface  at  k  =  L  —  k,T  +  1/2: 


Ak+1/2  ~  Ptop  +  (%+i/2  Bk+1/2)(pL+1/2  Ptop) 

B  k+1/2  ~  (bk+1/7)rk+1/2 


k  =  kp  +  1, . . . ,  L  -  ka. 


(A4) 


where  bk+V2  =  b(r\k+V2 )  is  defined  as  the  function: 


Hv) 


V  ~  Vkp+ 1/2 
1  —  Vkp+ 1/2 


(A5) 


and  the  exponent  r  in  the  Bk+ln  definition  in  (A3)  is  a 
function  of  fj:  rk+ll2  =  r( vk+1/2). 

Finally,  since  some  models  include  boundary  layer 


parameterizations  that  may  be  tuned  specifically  to  a 
levels,  the  algorithm  allows  the  lowest  ka  model  la¬ 
yers  to  be  cr-likc.  However,  use  of  the  pure  cr  coordinate 
[(4)]  for  these  levels  would  yield  a  discontinuity  in 
layer  thickness  at  the  L  —  ka  +  1/2  interface.  If 
rL_k  +i/2  =  1,  that  discontinuity  is  removed  by  setting 
these  ka  lowermost  levels  as  follows: 


Bu 


:  P  top 

/  Vk  +  l/2 


1  Pk  +  l/2  Bk+1/2 


)(Pl 


'  P  to 


Vkp+l/2 


\  1  Vkp  +  l/2 


\ k  =  L  +  1  —  ka,  ■■■,  L. 


(A6) 


These  surfaces  are  essentially  equivalent  to  a  cr- 
coordinate  that  levels  out  at  a  constant  top  pressure  of 
Pk  +1/2  instead  of  ptop:  more  formally,  they  are  pure 
Sangster-Arakawa-Lamb  (SAL)  levels  (see  section  3b). 
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